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ABSTRACT 

We have analyzed the motion of an infinitesimal mass in the restricted four body problem 
with solar wind drag. It is assumed that forces which govern the motion are mutual gravitational 
attractions of the primaries, radiation pressure force and solar wind drag. We have derived the 
equations of motion and find the Jacobi integral, zero velocity surfaces and particular solutions 
of the system. It is found that three collinear points are real when radiation factor < 3 < 0.1 
whereas only one real point obtained when 0.125 < (3 < 0.2. Again, stability property of the 
system is examined with the help of Poincare surface of section (PSS) and Lyapunov characteristic 
exponents (LCEs). It is found that in presence of drag forces LCE is negative for a specific initial 
condition, hence the corresponding trajectory is regular whereas regular islands in the PSS are 
expanded. 

Subject headings: Restricted four body Problem; radiation pressure; zero velocity surface; P-R drag;solar 
wind;PSS;LCE 



1. Introduction 

In space dynamics there are a number of sys- 
tems like two body, three body, four body, N-body 
problem etc. The simplicity and elusiveness of the 
three body problem in different form like restricted 
three body problem (RTBP), restricted four body 
problem (RFBP)(may be consider as an approxi- 
mation of two three body problem) etc. have at- 
tracted the attention of researchers for centuries. 
The motion of a spacecraft or satellite in the Sun- 
Earth-Moon system is a simple example of RFBP 
in space. The restricted four body has many pos- 
sible uses in the dynamical system for example, 
the fourth body is very useful for saving fuel and 
time in the traject ory transfers in the restricted 
four body problem ([Machuv et alll2007l) . 

The description of the effect o f radiation pres - 
sure force was first time given bv |Povntinj(| 19031) 
and the effect of total radiation force on a par- 



tic le P due t o radi ation source S was analyzed 
bv iRobertson (1937) with the help of the general 
relativity theory. He stated that if we consider 
only first order term in ^ then it consists a justifi- 
able approximatio n in classical mechanics to yield 
(jRagos et alill995h 
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(1) 



where F p — 16 ^^p Ac denotes the measure of the 
radiation pressure force, r is the position vector 
of P with respect to the Sun, v is the correspond- 
ing velocity vector and c is the velocity of light. 
In expression of F p , L is the luminosity of the ra- 
diating body whereas M, p and A are the mass, 
density and cross section of the particle respec- 
tively. In equation ([1} on right hand side, first 
term expresses the radiation pressure and second 
term represents the Doppler shift owing to the 
motion of the particle whereas third term comes 
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on account of the absorption and subsequent re- 
emission of the radiation. The last two terms are 
called P-R drag effect. 

Now, due to the sol ar wind drag force , equation 
(JT|) can be written as ( Burns et al 19791) 



mr 



SAQ f 



-(1 



1 - (1 + sw)- 



sw) 



(2) 



where S and Q pr are the solar energy flux den- 
sity and the radiation pressure coefficient respec- 
tively whereas sw represents the ratio of solar wind 
drag to Poynting-Robertson(P-R) drag (jLiou et al 
1995). The velocity independent term in ([2} de- 
notes the radiation pressure however velocity de- 
pendent term represents the drag force. The solar 
wind drag arises from the interaction between so- 
lar wind ions and the dust particles. The radiation 
pressure coefficient Q pr depends on the properties 
of the particle P and the radiation factor j3 is de- 
fined b 8 radiation pressure force SAQ pr r^ 

J ^ solar gravitat ion fo rce cGMra 

Few years ago, I Burns et al (1979) discussed 
about the radiation forces on small particles in the 
solar system and examined the differ ent types of 
effect of the radiating body. However. ISchuerman 
(1980) determined the equilibrium points and ex- 
amined their stability in the presence of radiation 
pressure and P-R drag forces. The dynamical ef- 
fect of general drag force (i.e. gas drag, nebu- 
lar drag, PR drag etc.) in the planar circular 
restricted thre e body problem was described by 
iMurravl (jl994l) . Also, he has examined the sta- 
bility of Lagrangian equilibr ium points using lin- 
ear approximation. Further, iLiou et all (|1995l) an- 
alyzed the effect of radiation pressure, P-R drag 
and solar wind drag on the motion of dust grains 
which is trapped in mean motion resonances with 
the Sun and the Jupiter in the restricted three 
body problem and foun d that all dust ga i ns or- 
bits are unstable. Again, Kalvouridis et al (2006) 
discussed the effect of radiation force due to pri- 
maries in the restricted four body problem using 
Radzievskii's model and noticed that there are 
some variations in the result which are unstable 
for all v alues of the parameters ass ume d by him. 
Furth er Jlshwar and Kushvahl (j2006l ) and lKushvahl 
(|2008l ). studied the restricted three body problem 
with P-R drag and examined the effect of P-R drag 



force on the motion of infinitesimal body. 

Lyapunov characteristic exponent (LCE) and 
Poincare surfaces of section (PSS) are efficient 
method to study the behavior of orbit around the 
equilibrium points. The LCE are very useful tool 
for the estimation of chaoticity of the trajectory in 
a dynamical system. Basically, it measure the rate 
of exponential divergence between neighborhood 
trajectories in the phase space. The ba sic con- 
cepts of LCE was given bv lOseledecl ( 1968 ) during 
the study of ergodic theory of dynamical system. 
The description of numeric al algorithm to calcu - 
late LCE was presen ted bv iBenettin et all ( 1980h 
and Froeschlel (|l984 ) . Whereas the determination 
of the stability regions of t he infinitesimal body 
was first time introduced bv IPoinca rc (1892) dur- 
ing the study of periodic orbit of the system. This 
is a good technique to study the nature of tra- 
jectory of infinitesimal body and also known as 
surface of se ction method. Further, this method 
was used by IWinterl (|2000l) to describe the loca- 



tion and stability of orbit in the restricted three 
body problem. 

A number of various type of four bod y mod- 
els h a ye been studied by many authors dHuane 



1981 



1968; Hadiidemctriou 1980; 



Michalodimitrakis 
Rosaevl 120111 : ICeccaroni and Biggsl 2012 ) 



etc. In these models, they have assumed either 
one or more bodies as radiating and minor body 
as very small in comparison to massive bodies. To 
preserve conservative character these models have 
been studied by ignoring dissipative forces like 
Poynting-Robertson(P-R) drag and solar wind 
drag. However, some of the research work in 
RFBP in addition with radiation pressure, P- 
R drag etc. have been performed by many au- 



thors (Kalvouridis et al 20061: Papadakis 200 



Machuv et all 120071 : iMedvedev and Perovl 



etc. But, little attention has been paid on the 
effect of solar wind drag such models. Therefore, 
we have considered a restricted four body problem 
with radiation pressure and solar wind drag. 

In order to know the nature of trajectory of in- 
finitesimal body in the proposed problem, we have 
determined first order LCE for the maximal evolu- 
tion tim e < t < 1000 with the help of algorithm 
given in lSkokosI (|2010l ). The PSS of the proposed 
problem is also obtained with the help of Event 
Locator Method. We have used Mathematica® 



(jWolframl 120031 ) for numerical and algebraic com- 
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putation. 

The formulation of the problem and derivation 
of equations of motion are presented in section ([2]) 
whereas in section @, we have obtained the zero 
velocity surfaces with different cases and location 
of the collinear equilibrium points is computed in 
section (TJ| . Section (|4.4p , we have discussed about 
the non-collinear equilibrium points however in 
section (|5|) we examine the stability behavior of 
the trajectory in the present dynamical system. 
Finally, we conclude the paper in section (J6j) . 




2. Formulation of the problem and equa- 
tions of motion 

We have considered a restricted four body prob- 
lem in which mj, i = 1, 2, 3 be the masses of three 
bodies such that 777,3 > mi > 777,2, and m be the 
mass of infinitesimal body. It is assumed that 7773 
be a radiating body with solar wind drag which re- 
volves around common center of mass of mi and 
7712 ■ The governing forces of the motion of in- 
finitesimal body are gravitational attractions of 
three massive bodies in addition with radiation 
pressure force, P-R drag and solar wind drag. It 
is also assumed that the effect of infinitesimal mass 
on the motion of the remaining system is negligi- 
ble. 

We use the canonical system of units, divid- 
ing all the distances by the distance between two 
primaries and dividing all the masses by the to- 
tal mass of the two primaries. Under these as- 
sumption, the mean motion of primaries is taken 
as unity. 

The masses and distances of the Earth, the 
Moon and the Sun are described as follows: mass 
of the Eaxth(M B ) = 5.98 x 10 24 fcg; mass of the 
Moon (Mm) = 7.35 x 10 22 kg; mass of the Sun(M s ) 
= 1.99 x 10 30 kg; distance between the Earth and 
the Moon is c?i = 3.844 x 10 5 fcm; distance between 
the Sun and the Earth is g?2 = 1-496 x I0 8 km. 
Again, the masses of the Earth, the Moon and the 
Sun in the canonical system are given as fj,E = 
0.987871, hm = *, M r*,_ = 0.012151 



M E 



Mm+Me 



Mm+Me 



and pL S = m ^ Me = 328900.48 respectively. 

Let he and fiM be the radii of the Moon and the 
Earth respectively and (x,y),(x E ,y E ),(x M ,yM) 
and {xs, ys) be the co-ordinates of the spacecraft, 
the Earth, the Moon and the Sun respectively (fig- 
ure [1]). 



Fig. 1. — Geometry of the problem 



Let x E = -fJ-M cos(t),y E = —fiM sin(i), ijf = 
H E cos(t),y M = n E sm(t),x s = R s cos(V>), ys = 
Rs s'm(ip) and ip = ip + uj s t, where Rs = 
389.1724 is the distance between the Sun and 
the center of mass of the system, ujs — 0.07480 
is the angular velocity of the Sun which makes 
an angle ip with x-axis, ipo is the initial value 
of ip and t is the time. The distances of 
spacecraft from the Earth, the Moon and the 
Sun are given as r\ = y/{x — xe) 2 + (y — Ve)" 2 
t% = y/Qc - x M ) 2 + (y - ynY and r 3 = 
y/(x - x s ) 2 + {y- ys) 2 respectively. 

The equations of motion of the infinitesimal 
mass in the inertial reference system are 



where 



Fpr,x - 
FpR, v 



He(x - x E ) _ fHfjx - x M ) 

rf r\ 
/js(1 - (3)(x - x E ) 
r| 

+(l + sw)F PRlX , 
v<E{y - ve) _ hm{d - vm) 

r\ r\ 

Ms(i-/3)(y-ya) 
+(l + sw)F PR>y , 

Pv-s 



(3) 



x — xs) + X , 



-^r [h(y - ys) + y] 

cri 



(4) 

(5) 
(0) 
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But, 



{(x - x s )x + (y- ys)y} 



In the above expression qi = (1—0) — (1— -jr) and 
c are the mass reduction factor and dimensionless 
speed of the light respectively whereas F p and F g 
are the radiation pressure and the gravitational 
force respectively. 

Using the coordinate transformation x = 
£ cos t — rj sin t; y — £ sin t + r\ cos t, equations ^ 
and ^ become in the rotating reference frame as: 



£ - 277 = £ - 



Mi-/9)(£-&) 



+ (1 + aw)Fpn,e, (7) 



V + ^^r, -3 -3 



or 



at/ 



(1 + fliuJFpfl,,, (8) 



t - 2v? = -jjr + (1 + sw)F PRti , 



dU 

V + ^= -7^ + (l + sw)F PRiV , 



where 



jj _ Re _|_ /Aa/ _|_ /Jg(l-/3) + £ 2 + t; 2 



n r 2 



''3 



(9) 



(10) 



(11) 



is the potential in the rotating coordinate system. 
The P-R drag components in (£, rj) reference sys- 
tem are 



Fpr,£ - 
and 
FpR,n 
where 



fas 

2 



[i - (V - Vs)} + ni(£ - &) 



[{v + (£- Cs)} +ni(rj - ?7 S )] , 



"i 



{(e-^)£ + (r?-r7s)r)} 



If we multiply the first equation © by £ and the 
second equation (fT0|) by t) and adding, we get 



(1 + sro) 



9E7 dS, + dU dr) 
d£ dt drj dt 
~dFd£ + dF drj 
dt; dt drj dt 



(12) 



Since U and F do not depend explicitly on time, 
the expression on the right hand side is the total 
time derivative of U and F. The left hand side 
can be expressed in terms of the derivatives of the 
velocities: 



1 d r 

2 ~dt 



dU 
~dt 



(1 



SKI 



dF 
It 1 



(13) 



where 

dF 
~dt ' 



fas 



{(£-Zs)i + (v-vs)fi} A 



-{(v-Vs)£-(Z-Zs)ri}+£ 2 +V 2 ] ■ (14) 

Integrating of equation (fl~3f with respect to time, 
gives 

£ 2 + r) 2 = 2U + 2(1 + sw)F - C 



or 



C = 2U + 2(1 + sw)F 



(15) 



where v is the velocity of infinitesimal body and 
C is a constant, called the Jacobi integral. Now, 
the square of the velocity cannot be negative, the 
motion of the negligible body is restricted to the 
region where 



v 2 = 2Z7 + 2(1 



or 



U 



(1 + sw)F > — 



)F-C>0 



C 



(16) 



When the position and velocity of the negligible 
body are known for some initial moment, then 
we can obtain the Jacobi integral. Since U and 
F are functions of position and potential of dissi- 
pative force only, condition (|16|) tells immediately 
whether the system can ever reach a given point (£, 
77). This condition does not tell about the shape 
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of the orbit, it only determines the region where 
the particle could move. Also, condition (fTB| shows 
that the larger value of C, the smaller the allowed 
region. 

When C is large, the allowed region consists 
of the four separate areas. It can never moves 
form one allowed region to another if the regions 
are not connected. When C becomes smaller, a 
connection opens first at the point L\. Again, we 
take C even smaller, a connection opens between 
two inner regions, first at the point L2 and second 
at the point L3. The small body can never escape 
from the system. Its orbit is then stable. Finally, 
when the value of C is further reduced, the outer 
and inner regions become connected and escape 
is possible which can be seen in the figure (|2|). 
It is also observe that how does the connection 
open with the decreases value of C . Figure© and 
(jj) are the zoom portions of the regions around 
the Earth and the Moon. The points L\, L2, L3 
lie on the same straight line which connects the 
primaries. 

When there are no dissipative forces i.e. Fprj = 
0, then constant of the motion and the Jacobi con- 
stant C is defined as 



C = 2U-U 2 +f, 2 ), 



(17) 



which gives the results as in classical model of re- 
stricted three body problem. Again, equation (|15|) 
can be written as 



where 



C = 2U + h + b 2 - - if 



(18) 



1 



+arctan 
62 = 



r l2((£- 389.172) 2 + 77 2 ) 

n 



£- 389.172 

2#*sr(l + ow) 



and 

(i 2 + v 2 ) 



- 389.172) 2 + ?7 2 



dt. 



The third term of equation (|T8|) depends on the 
time due to drag force. Therefore this integration 
term shows that Jacobi constant depends upon th e 



time ([Murray and Dermottlll999HLiou et ajll995t ). 

Consequently, with the help of equation (fTT|) . 
we obtain the time derivative of Jacobi constant 
in presence of drag forces which is given as 



C = 



2(1 + sw) — 2" 



{{ri-Vs)i-(£-ts)v}} (19) 
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Fig. 4. — Zoom part of figure © at £ 
(0.9, f.f), r\ <= (-0.1,0.1), Ci = 848.255, C 2 
742.489, C 3 = 679.026. 



aries of regions in which the particle is free to 
move. Now, for large values of £ and r], all the 
terms except first and second in L.H.S. of the equa- 
tion (|20p become unimportant. In other words, 
this equation takes the form: 



e+rf =C-a = d, 



(22) 



where 



2a M , 2(1 - P)fi s 

1 h A\. 



r-2 



T3 



The equation (|22l) represents a circle whose radius 
is \/Ci. Therefore the curve in the £77- plane is 
an approximately oval shape within the asymp- 
totic cylinder. For small values of £ and r], the 
first and second terms are relatively unimportant, 
hence the equation may be written as: 



Mm (1 - P)ns 
r-i r 3 



A! = C 2 



where 



We have tried our best to derive the explicitly time 
independent expression of potential function or Ja- 
cobi constant but we could not find the exact ana- 
lytical expression. In future, this should be taken 
up as a challenging problem of such models having 
P-R drag and solar wind drag. 

3. Zero velocity surfaces 

The Jacobi integral is a relation between the 
square of the velocity and the coordinates of the 
infinitesimal particle with respect to the set of ro- 
tating axes. If the particle's velocity becomes zero, 
we have 



e+v 2 

where 
A, = 



2^_£ 



2/j M 2(1- P)n s 



ri 



2(1 + sw)pn s 



'"3 



arctan 



A x = C, (20) 



£-£s 
V~Vs 



-arctan 



( V 


- nsX 


u 





(21) 



and C can be determined from the initial condi- 
tions. 

Equation (|20[) is important in this problem and 
it is defined for a given value of C i.e. the bound- 



(23) 



which is an equation of equipotential surfaces. 

Figure ([5]), shows that curves intersect at a 
point where function gives local minimum and 
maximum value of Jacobi constant. Again, figure 
(J6j> shows contour curves which are the equipoten- 
tial boundaries. It is clear from the frames that the 
configuration rotates anti-clock wise with an angle 
ifi. In the contour plots red (dark) portion shows 
the higher equipotential value whereas when we go 
from red portion to blue (very dark) one, equipo- 
tential decreases gradually. Now, the zero velocity 
curves are defined by the equation 



2fi E 2/j m 



+ 



n r 2 

2(1 



''3 



(24) 



To determine the curvature of a zero velocity 
curve (ZVC) (|2~4"1) , we use the formula (|Marceta 
2012b 



zvc ~ UfTJIW 2 ' (25) 



6 



Fig. 6. — Zero velocity surfaces at (a) ip = 0, (b) 
ip = 60, (c) tp = 90 and (d) t/> = 180. 



Fig. 5.— (a) f(£, 0) and (b) f(0, 77) for /i E 
0.987, fi M = 0.012, ^ s = 328900.48, t 
1000, V = 0, P = 0.00002, c = 1.8 x 10 6 , sw 
0.248. 



where 



, df^rj) dfforj) d 2 f(^ V ) 



drj 1 



hi 



dr) 

d 2 f(t;,v) 
aft • 



Now, we determine the curvature of ZVC (124)) 
with the help of characteristic values of Hessian 
matrix H, i.e. 



H 



dN ( dN ( 
Of dn 

d£, dr] 



(26) 



where N% and N v are components of normal unit 



vector defined by 



N = 



Of df_ 



G_ 

W\ 



1/2 : 



(27) 



with 



G = v/ = 



a/ a/ 



(28) 



Since the expression of the curvature of zero ve- 
locity curve is too complicated to deduce anything 
meaningful analytically. Therefore we have com- 
puted curvatures numerically and observe that the 
curvature of any curve shows the bending nature 
at a particular point (figures [7] , [5] and H]). This 
value is inverse proportional to the radius of cur- 
vature. If curvature is large then radius of cur- 
vature is small, that means curve is more bended 
at that point. If curvature decreases than radius 
of curvature increases, that means bending prop- 
erty of the curve decreases which opens the path 
to move particle freely. From figure© it can be 
seen that around the singular points, curvature of 
zero velocity curve is highest that means radius 
of curvature is least therefore curves are disjoint 
closed loops form where particles can move form 




Fig. 7.— Curvature of (a) f(f , 0) and (b) f(Q, rf) 

for (is = 0.987, /i M = 0.012, ^ s = 328900.48, t = „. . _ , , . . 

0, V = 0, p = 0.00002, c = 1.8 x 10 6 , sw = 0.248. Flg ' 9 ~ Curvature of zer0 velocit y surfaces - 



one regions to another. In figure ([7]). frames (a) 
and (b) show the curvature of ZVC depends on 
£ and ry respectively and indicates that the curves 
are sometimes nearer to each other and sometimes 
goes far form to each other. When the curves bend 
to each other then the connection opens where the 
body moves from one allowed side to another side. 
When the curve is not bended that means curve 
goes outside, then the curves are separate (figure 
[5]). Different layers show the equipotential cur- 
vatures where the movement is possible which is 
clearly seen in zoom part of the figure jSJ). Also, 
figure shows the curvature of zero velocity sur- 
faces, cut off portions indicate the singularities of 
the surfaces where motion of infinitesimal body is 
impossible. 



4. Location of the Lagrangian equilibrium 
points 

For equilibrium points, we solve the equations 
© and (flOj by taking £=0 = r? = £ = ?7 i.e. 

t _ Mg(£ - gg) _ mm(£-£m) 

? r 3 r 3 

'1 '2 

3 h (1 + sw)F PRt £ = 0, (29) 

Me ~ Ve) Mm (v - Vm ) 

7 1 3 3 

ms(i - P)(v - Vs) , ^ , n , qn s 
3 h (1 + sw)bp R<n = 0, (30) 

r 3 

where 

F PR<i = eg-iv-Vs), F PRtn = 



s 



4.1. Particular cases 

1. When = then the equation ((29j and (|30|) 
give the following 

> _ ve{Z-£,e) _ MHm) _ ms(£ - £g) = n 

5 _a _a _a > 



He{v~Ve) (J-m(v-Vm) Vs{v~Vs) 



Now, if we take j3 non-zero then we obtain three 
real roots of equation (|33|) change with the value 
of j8 G (0,0.1). Again, if £ e (0.125,0.2) then we 
have only one root which increases with (3. With 
the help of these values we determine Jacobi con- 
stant and see that an increases in the value of 
/3 consequently the Jacobi constant decreases and 



This is the classical case of restricted four body 
problem. 

2. If = 1 then the equation and ^ 
takes the form 



= 0. (31) the corresponding energy increases (Tabled]) 

4.3. Extremal values 

From equation (|20j) . we have 

2/lE 2jU M 



V - 



he^-Ze) mm(£-£m) 



\ie{v - Ve) Mm (v - Vm ) 



+ (1 + sw)F PRX = 0, 
+ (1 + 8w)F PR>7l = 0. (32^, 



C = £ 2 + TJ 

2(1 + sw)^ s 
c 

ow, 



(»7 - Vs) 



r-2 



2(1 



(35) 



This becomes the restricted three body problem 
with solar wind and P-R drag. 

4.2. Collinear points 

The collinear points are solutions of the equa- 
tions (|29H30p with -q = 0. That is, the collinear 
points lie on the line joining the primaries. There- 
fore, we have 



He((,-£.e) Mm(£-£m) Ms(f-£s) 



0. (33) 



When the value of j3 is zero, the above equation 
becomes 



£ - 780. 294£ 6 + 152974£ 5 
-624960C 4 + 782507£ 3 - 454193£ 2 
+288215^- 145671 = 0, 



(34) 



which shows that the roots lie in the £ axis and 
after simplification, the above equation consists 
seven degree. Using the Descartes' rule of sign we 
found that three roots are real and other pair of 
roots are imaginary, because the coefficients of the 
power of £ change seven signs. This shows that at 
a time we can have only three real roots joining the 
two primaries because the plane of motion of the 
Earth around the Sun is different from the plane of 
motion of the Moon around the Earth. The three 
primaries will be in the same plane when the Moon 
comes at the line of node of the plane of motion. 



HevI + fiMrl = he [(£ - £e) 2 + (v- Ve) 2 ] 
+Mm [(£ - Cm) 2 + (v- vm) 2 ] , 
=> m v \ + mm^I = £ 2 + v 2 + he^m, 

=> £ 2 + T] 2 = HET\ + l J, M r 2 - VEV-M- 

Substituting these values in equation (|35[) . we get 



Mb 



2 
ri 

2(1 + aw)Pns 



Mm 



2 

''2 



J(v-Vs)d£ f(t-£s)dv 



-MsMm 



2(1 -/%s 



T3 



(36) 



Differentiating equation (I36|) with respect to 7*1 , r 2 
and r% respectively and equating each term to zero 



Table 1: Collinear 
He = 0.987, fi M = 
ip = 



points for initial parameters 
0.012, ^ 5 = 328900.48, t = 0, 








I 


II 


III 


Cl 


c 2 


C3 




0.0 


-1.86495 


-0.892772 


0.578692 


843.374 


844.720 


848.255 





0000025 


-1.86494 


-0.892775 


0.578693 


843.372 


844.718 


848.253 


0.000005 


-1.86494 


-0.892778 


0.578693 


843.370 


844.716 


848.251 





0000075 


-1.86493 


-0.892780 


0.578694 


843.368 


844.714 


848.249 




0.00001 


-1.86492 


-0.892783 


0.578694 


843.366 


844.712 


848.247 





0000125 


-1.86491 


-0.892786 


0.578695 


843.364 


844.709 


848.245 





0000150 


-1.86490 


-0.892789 


0.578695 


843.362 


844.707 


848.242 





0000175 


-1.86490 


-0.892792 


0.578696 


843.359 


844.705 


848.240 




0.00002 


-1.86489 


-0.892795 


0.578696 


843.357 


844.703 


848.238 




0.025 


-1.78593 


-0.924049 


0.583721 


822. .399 


823.571 


827.103 




0.05 


-1.70064 


-0.961693 


0.588858 


801.429 


802.416 


805.950 




0.075 


-1.60522 


-1.009590 


0.594104 


780.467 


781.248 


784.796 




0.10 


-1.48925 


-1.078140 


0.599462 


759.527 


760.056 


763.643 




0.125 


Complex 


Complex 


0.604934 


Complex 


Complex 


742.489 




0.150 


Complex 


Complex 


0.610524 


Complex 


Complex 


721.335 




0.175 


Complex 


Complex 


0.616232 


Complex 


Complex 


700.181 




0.20 


Complex 


Complex 


0.622061 


Complex 


Complex 


679.026 
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for initial values of t, we get 



ri 



1, r 2 = 1, 

/3 2 x (1 + sw) x R s x ry) 
/3 - 1 c ' 



(37) 



In equation (|37j). if we take r3 = 0, this causes 
the Sun to be placed at the origin which is not 
possible in the proposed system. Because for a 
moment if we consider that the Sun is placed at 
the common center of the Earth and the Moon and 
consequently, we get a system in which the Earth 
and the Moon revolve in different orbits about the 
Sun. In general such type of system does not exist 
for the problem in hand. If r$ > then the factor 



> has the following cases: 



(l+sw)P 

1. If rgZ/n (1 + sw ) have positive signs then 
sw > — 1 and /3 lies in the open interval (0, 1) 
that is |? < 1. 

2- If (/3-i) ana - (1 + sw ) have negative signs 
then sw < — 1 and /3 lies outside the open interval 
(0,1) and the body escapes from the system. 

The function </> always positive for any rj, i = 
1,2,3 and </> approaches to infinity as r, approaches 
to zero or infinity. Therefore, the equilibrium 
points of the system are absolute minimum how- 
ever the pseudo-potential and the Jacobi constant 
at these points are: 



1 



(1 + sw) /3 rj 



-4, 



C 



(1 + sw) (3 rj 



-An 



(38) 



(39) 



where 



A 2 = -7.60616 x 10 s + j3 {1.52123 x 10 9 

+ (2.98516 + 2.98516 x sw)} . (40) 

4.4. Non-collinear points 

For non-collinear points, we take perturbation 
in distances of infinitesimal mass from the pri- 
maries due to the radiation pressure and solar 
wind drag and suppose rj = 1 + ei and r 2 = l + e 2 , 
where ei,C2 <C 1, then from center of mass prop- 
erty, we obtain 

V-E-ri + ^M r 2 = =» He(\ + A*M e 2 + 1 = 0. (41) 



Now, r\~rl = (£+Mm) 2 +V 2 - (£ + Mm - I) 2 -?? 2 , 
therefore 



2(ei-e 2 ) + l 
? = ^ Mm- 



(42) 



Substituting this value in the equation <f29|) . we 
get 

/ -(l-3aw-4/3) xc 2 x (g - 389.172) 
V V 2 x (1 - sw - 2/3) x 389.172 x p 2 ' [ ' 

Again, substituting the value of £ in the equation 
([29)) and neglecting second and higher order terms 
of /3, we get 



/i S (l-4/3)(l-3su;)c : 



3 , 2(ei - e a ) + 1 



-Mm - 389.172) = 0. (44) 



Solving equation (|4Tj) and (|44| . we get the values 
of ei and e 2 which are given below 



1 



ei = 



£2 = 



-l + 3sio + 4/3 d ' 



1 



-l + 3sw + 4/3 



tA 4 



(45) 



where 



A, 



and 



-12(0.0833 - 0.3333,3 - 0.25sw) x 
(0.0026 + fi M ) x (388.669 + fj, M ) (46) 



A± = 1 - 3sw - 4/3 + (4664.06+ 12^ M )x 

^(0.0833- 0.3333/3 - 0.25sw). (47) 

Substituting the values of e\ and e 2 in equa- 
tion (|4"2"|) and (|4"3"|) , we get the non-collinear points. 
These points depend on the values of the radi- 
ation pressure and the solar wind drag. When 
we change the values of ratio of solar wind to 
P-R drag parameter for different fixed values of 
radiation pressure parameter, points are changes. 
Numerically the coordinates of triangular equilib- 
rium points are L 4 f = (0.615824, 0.75845), L 5 f = 
(0.615824, -0.75845) which are obtained by tak- 
ing the same parametric values as earlier taken in 
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-60000 -40000 -20000 20000 40000 60000 
fit] 

Fig. 10. — Trajectory of non-collinear points 



case of collinear equilibrium points. The paramet- 
ric solution (i.e. triangular point) is depicted in 
figure (fTU|) . Frame (a) of this hgure shows that 
triangular point exists regularly for < t < 60 
and then a sudden change in trajectory has been 
seen for 60 < t < 80. Finally, solution maintains 
its regular path for t > 80. Whereas frame (6) 
of this figure shows the nature of the path of in- 
finitesimal body for the time interval < t < 200. 
It is also found that the infinitesimal mass escapes 
out side of possible region after t > 424229. 

5. Poincare surfaces of section and Lya- 
punov characteristic exponents 

The Poincare surfaces of section (PSS) is an 
important technique to describe the stability re- 
gion of the system. In order to determine the or- 
bital element of the infinitesimal body at any in- 
stant, it is necessary to know its position (£, 77) 
and velocity (£, fj) correspond to a four dimen- 
sional phase space. In this paper, we have deter- 
mined PSS in the ££-plane using Eyent L ocator 
Method of Mathematica® (|Wolfram| |2003|) . The 
magnitude of the velocity component f] is com- 
puted with the help of Jacobi integral equation 
(fT5|) however the Jacobi constant C initially com- 
puted at t = 0. Now, we have plotted the graph 
of £ and £ against each other when the trajectory 
intersects the plane in the direction of 77 > 0. In 
other words, a smooth well defined island in the 
PSS indicates that the trajectory is likely to be 
regular whereas the fuzzy distribution of intersec- 
tion points represents chaotic trajectory. Again, 
if the curve shrink to a point that means it has 
a periodic orbit. Further, we have obtained PSS 
at the values of Jacobi constant C for a certain 
range of values of £ and £ whereas each orbit is 
determined with initial conditions: 

( = Zo, r) = Q, i = 0, 



r)= ^h + b i + e a ~e Q -C, (48) 



where 

63 = 
and 



2/-iM 



2Ml-j8) 

(Co + 0.012) ^ (Co - 0.987) (Co - 389.172) 



In 



1 



-arctan 



L2(£o- 389.172) 2 

n 



Co - 389.172 
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Since, in the proposed model key quantities are the 
values of C, f3 and sw. Therefore, we have plotted 
the graphs at these parameters. Figure (JTTJ) shows 




it is clear that the regular island expands gradually 
due to the effect of radiation pressure and solar 
wind drag. Again, the island at £ = 0.025 shows 
that the trajectory is regular i.e. region in the 
neighborhood of £ = 0.025 is stable. On the other 
hand, in frame (b) the island at 0.025 shows that 
the trajectory is regular and this region shrinks 
towards center. Hence from frames (a) and (b), 
we conclude that the radiation pressure and solar 
wind drag have significantly affect on the stability 
region of the trajectories. 



- -1.0 

B 

J -1.5 
-2.0 
-2.5 



(a) 




Fig. 12. — Lyapunov characteristic exponent: (a) 
for < t < 10000 (b) Zoom of (a). 



Fig. 11. — Poincare surfaces of section. 

PSS for the restricted four body problem at the 
minimum values of Jacobi constant C — 1696.48 
which is initially determined at £ = 0. 5,^ = 0,^ = 
0, fj = 0. In this figure, we observe that the nature 
of PSS and the trajectory in ££-plane. From frame 
(a), we have seen regular island with the radiation 
pressure and the solar wind drag but frame (b) 
shows the regular island without any force. Here, 



In order to know the behavior of nearby trajec- 
tory emanating from the neighborhood of an equi- 
librium point we compute the LCEs. With the 



help o f method and algorithm described in lSkokos 
(2010), we have numerically computed first order 
LCEs and plotted the graphs t Vs LCE(A), upto 
maximum evolution of time < t < 10000 when 
He = 0.987, (j, M = 0.012,^s = 328900.48,/? = 
0.00002, sw = 0.248, c= 1.8xl0 6 (figurc[l2]). This 
figure shows that the proposed model is slightly af- 
fected in presence of the radiation pressure and the 
solar wind drag. Also, it is clear that A decreases 
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slowly with time. Furthermore, A is negative for 
< t < 10000 which shows that orbit is regular 
in this interval. Therefore, we observed that the 
values of LCE depends on initial deviation vec- 
tor and renormalized time step. As we have used 
default method of Mathematica® software pack- 
age for the numerical solutions of the differential 
equations. The obtained result depends on preci- 
sion and accuracy goals, i.e. in present problem, 
if these goals were chosen different from 3 con- 
sequently there was overflow in the result which 
leads to the false estimation of the LCE. 

6. Discussion and conclusion 

We have studied the restricted four body prob- 
lem by assuming the effect of radiation pressure 
and solar wind drag. The boundaries of allowed 
regions for the motions of the infinitesimal mass 
are determined using zero velocity surfaces at dif- 
ferent values of radiation pressure and solar wind 
drag. It is found that allowed possible regions of 
the motions decrease with the an increase in the 
value of Jacobi constant C. We have observed that 
in presence o f drag forces, Jacobi constant de pends 
on the time (|Murravlll99l iLiou et allll995l ). The 
range of radiation pressure and solar wind drag 
are determined and found that motion is possi- 
ble for the values lie in the interval < j3 < 1 
and — 1 < sw < 1 respectively. The curvature of 
the ZVC is obtained and noticed that when the 
curves bend to each other then connection is open 
and the body moves from one allowed side to other 
side however when the curve is not bended then 
curve goes outside and the body does not move 
one allowed side to other side. 

We have obtained the particular solutions, 
which are the extreme values of Jacobi function. 
It is found that three collinear points for the val- 
ues of < (3 < 0.1 are real whereas only one real 
point exists when 0.125 < (3 < 1. We have deter- 
mined the co-ordinates of two non-collinear points 
which depend on radiation pressure and the solar 
wind drag. With the help of PSS, it is observed 
that the stability region get expanded in presence 
of radiation pressure and solar wind drag and at 
the point £ = 0.025 orbits are stable. Further, 
we estimated the LCE(A) for the maximal time of 
evolution < t < 10000 and found that it is neg- 
ative which shows that orbit of the infinitesimal 



body is regular. Since it is difficult to obtained an 
exact expression of pseudo potential function in 
presence of term due to dissipative force, therefore 
further work is needed in this regard. This work 
may be applicable to study the motion of test 
particle in coupled restricted three body problem 
with drag forces. 
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